##################################################################
# function normalize
##################################################################w
normalize <-function(matrix) {
	# dimensions of the matrix
	n1<-dim(matrix)[1];
	n2<-dim(matrix)[2];

	#pre-normalize in the 1-direction
	ssquare<-apply(matrix*matrix,1,sum,na.rm=TRUE);
	#repmat: repeat and reshape
	sumTerms<-rep(ssquare,n2);
	dim(sumTerms)=c(n1,n2);
	normalize<-matrix/sqrt(sumTerms);	

	#subtract mean in the 2-direction
	sumTerms<-apply(normalize,2,mean,na.rm=TRUE);
	#something like repmat: repeat, reshape and transpose
	sumTerms<-rep(sumTerms,n1);
	dim(sumTerms)=c(n2,n1);
	sumTerms<-t(sumTerms);
	normalize<-normalize-sumTerms;
	
	#set sum(squares) = 1 in the 2-direction
	ssquare<-apply(normalize*normalize,2,sum,na.rm=TRUE);
	sumTerms<-rep(ssquare,n1);
	dim(sumTerms)=c(n2,n1);
	sumTerms<-t(sumTerms);
	normalize<-normalize/sqrt(sumTerms);

	# replace missing values by zero
	normalize[is.nan(normalize)]<-0;
	return(normalize);
}


##################################################################
# signature algorithm simple application
##################################################################w
sigAlg<-function(genesIn,nd1,nd2,threshC,threshG) {

	#number of genes in the input set
	nGenesIn<-apply(genesIn,2,sum);
	
	#number of input sets
	nInputSets<-dim(genesIn)[2];

	#data dimensions
	nGenesGenome<-dim(nd1)[1];
	nConditions<-dim(nd1)[2];

	#calculate expected standard deviation
	thresholdConditions<-threshC/sqrt(nGenesIn*nGenesGenome);

   	#-------------------------------------------------------
	# First step of algorithm: Score conditions
	#-------------------------------------------------------
	conditionsOut<-t(genesIn)%*%nd1;

	#divide by number of genes to get average
	nG<-rep(nGenesIn,nConditions);
	dim(nG)=c(nInputSets,nConditions);
	conditionsOut<-conditionsOut/nG;
	
	# apply threshold
	meanMatrix<-rep(apply(conditionsOut,1,mean),nConditions);
	dim(meanMatrix)=c(nInputSets,nConditions);
	thresholdMatrix<-rep(thresholdConditions,nConditions);
	dim(thresholdMatrix)=c(nInputSets,nConditions);
	
	conditionsOut[abs(conditionsOut-meanMatrix)<thresholdMatrix]<-0;

   	#-------------------------------------------------------
	# Second step of algorithm: Score genes
	#-------------------------------------------------------
	genesOut=nd2%*%t(conditionsOut);
		
	# apply threshold	
	# for genes, work with MEASURED standard deviation
	# which is essentially same as the expected
	meanMatrix<-rep(apply(genesOut,2,mean),nGenesGenome);
	dim(meanMatrix)=c(nInputSets,nGenesGenome);
	meanMatrix<-t(meanMatrix);
	thresholdMatrix<-rep(threshG*apply(genesOut,2,sd),nGenesGenome);
	dim(thresholdMatrix)=c(nInputSets,nGenesGenome);
	thresholdMatrix<-t(thresholdMatrix);
	
	genesOut[(genesOut-meanMatrix)<thresholdMatrix]<-0;

	return(list(genesOut,t(conditionsOut)))
		
}



